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Abstract 

We  report  the  results  of  extensive  numerical  experiments  on  the  Kuramoto-Sivashinsky  equation 
in  the  strongly  chaotic  regime  as  the  viscosity  parameter  is  decreased  and  increasingly  more  linearly 
unstable  modes  enter  the  dynamics.  General  initial  conditions  are  used  and  evolving  states  do  not 
assume  odd-parity.  A  large  number  of  numerical  experiments  are  employed  in  order  to  obtain 
quantitative  characteristics  of  the  dynamics.  We  report  on  different  routes  to  chaos  and  provide 
numerical  evidence  and  construction  of  strange  attractors  with  self-similar  characteristics.  As  the 
“viscosity”  parameter  decreases  the  dynamics  becomes  increasingly  more  complicated  and  chaotic. 
In  particular  it  is  found  that  regular  behavior  in  the  form  of  steady  state  or  steady  state  traveling 
waves  is  supported  amidst  the  time-dependent  and  irregular  motions.  We  show  that  multimodal 
steady  states  emerge  and  are  supported  on  decreasing  windows  in  parameter  space.  In  addition  we 
invoke  a  self- similarity  property  of  the  equation,  to  show  that  these  profiles  are  obtainable  from 
global  fixed  point  attractors  of  the  Kuramoto-Sivashinsky  equation  at  much  larger  values  of  the 
viscosity. 
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1  Introduction 


The  Kuramoto-Sivashinsky  equation  is  one  of  the  simplest  one-dimensional  PDE’s  which  exhibits 
complex  dynamical  behavior.  As  an  evolution  equation  it  arises  in  a  number  of  applications  includ¬ 
ing  concentration  waves  and  plasma  physics  ([4],  [20],  [21],  [22]),  flame  propagation  and  reaction 
diffusion  combustion  dynamics  ([-30],  [31]),  free  surface  film-flows  ([2],  [14],  [29])  and  two-phase 
flows  in  cylindrical  or  plane  geometries  [7],  [11],  [25],  [35],  [36].  The  equation  can  be  written  as 


Ut  +  UUx  +  Uxx  "b  Ux 


*  ^  L  o  ’  o  ’ 


noting  that  the  solution  remains  periodic  of  period  L  if  the  initial  data  are  periodic. 

Equation  (1)  can  be  normalized  to  27r-periodic  domains  by  the  transformations  xi  =  ■Kxj L, 
U  =  iiu/L  and  ti  =  which  on  dropping  the  subscripts  1  yields  the  scaled  KS  equation 


Vtf  -|-  ttUx  T  '^xx  "b  ^y^xxxx  — 

{x,t)  6  IR^  X  IR"'’, 
u(x,  t)  =  u(x  +  2ir,  t). 


(2) 


where  i/  =  ^  >  0  represents  a  dimensionless  wavelength  of  u.  Physically,  then,  the  limit  — v  0 
corresponds  to  the  waves  in  (2)  becoming  infinitely  long.  It  is  easily  established  by  linearization  of 
(2)  that  for  a  given  i/  the  first  mod(i/~2)  Fourier  modes  are  linearly  unstable  and  grow  exponen¬ 
tially.  It  is  well-known,  from  numerical  experiments  for  instance  (see  below),  that  solutions  to  (2) 
become  increasingly  irregular  as  u  is  decreased. 

Without  loss  of  generality  the  study  of  solutions  with  zero  spatial  mean  is  sufficient  due  to 
Galilean  invariance.  An  analytical  study  of  the  KS  equation  was  carried  out  by  Nicolaenko  et  al. 
[24]  (referred  to  as  [NST]),  where  it  is  shown  that  if  the  initial  data  are  in  and  are  of  odd-parity 
(i.e.  antisymmetric  about  a;  =  0  for  (1)  or  x  =  ir  for  (2)),  then  the  solutions  remain  in  for  aU 
time  and  there  is  a  globally  attracting  set  also  bounded  in  L^.  The  bound  depends  on  the  value  of 
L  or  V  and  the  estimate  of  [24]  gives  the  following  upper  bound: 


lim  sup||u(-,t)||2  <  const  •  12  (3) 


The  corresponding  bound  for  (1)  is  0(i®/^).  Such  bounds  of  the  norm  are  useful  in  determining 
estimates  of  the  Hausdorff  dimension,  dji,  of  the  universal  attractor.  In  fact  working  with  equation 
(1),  it  has  been  shown  by  Temam  [34],  that  if  the  bound  for  the  norm  is  0(T“)  then  djj  is 
proportional  to  and  so  the  [NST]  result  yields  the  estimate  for  d}j;  we  note  that 

the  conjectured  best  bound  is  0(L)  but  a  proof  is  not  available  yet. 

The  analysis  in  [NST]  was  based  on  the  odd-parity  of  the  initial  data.  This  restriction  has 
been  removed  recently  by  several  authors.  Goodman  [13]  considers  general  smooth  initial  data  and 
obtains  a  bound  of  the  same  form  as  (3).  This  bound  is  improved  further,  however,  by  H’y^shenko 
[17]  and  independently  using  a  more  classical  approach  by  Collet  et  al.  [5],  and  is 

lim  sup||u(-,t)||2  <  const  •  (4) 

t—yco 


Again  the  analogous  bound  for  (1)  is  improving  the  estimate  of  the  Hausdorff  dimension 

to  bringing  it  closer  to  the  conjectured  best  bound.  We  note  that  the  analysis  in  [5]  and  [17] 

applies  also  to  odd~parity  initial  conditions  and  so  is  an  improvement  of  tie  method  in  [NST]. 

Analyticity  of  the  solution  has  also  been  proved  by  Collet  et  al.  [6].  In  fact,  considering  equation 
(1),  the  result  proved  is  that  for  large  the  function  U{x^t)  is  analytic  in  ic  in  a  strip  of  width 
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13  >  const  •  around  the  real  axis,  which  in  turn  implies  that  the  high  frequency  part  of  the 

spectrum  hajs  the  form 

\Unit)\  Ri  O(exp(-const  •  (5) 


where  Unit)  is  the  nth  Fourier  coefficient  of  U{x,t)  and  q  =  27r/i.  By  a  series  of  numerical 
experiments  a  much  stronger  result  that  (5)  has  been  conjectured  in  [6],  namely  that  there  exists 
a  /?  >  0  independent  of  L  such  that  the  solutions  of  (1)  satisfy 


lim  sup  Y]  ^  <  const., 

t-i-OO 


(6) 


and  the  numerical  work  gives  /?  «  .3.50.  These  results  are  extremely  useful  in  guiding  our  numerical 
work  as  is  explained  in  Section  2. 

There  are  several  computational  and  analytical  studies  dealing  with  the  construction  of  approx¬ 
imate  inertial  manifolds  of  KS.  Kirby  [19]  develops  a  Galerkin  approximation  based  on  Sobolev 
eigenfunctions;  Chen  [3]  constructs  approximate  inertial  manifolds  by  broadening  the  gaps  in  the 
spectrum  in  order  to  obtain  low  dimensional  behavior;  Foias  et  al.  [10]  construct  fully  discrete 
nonlinear  Galerkin  schemes  based  on  approximate  inertial  manifolds  of  KS;  Robinson  [28]  uses  the 
bounds  of  [5]  and  [17]  to  produce  a  new  estimate  of  the  dimension  of  the  globally  absorbing  set  of 
KS.  Additional  references  may  be  found  in  these  articles. 

A  significant  amount  of  numerical  work  preceded  the  theoretical  results  outbned  above.  The 
numerical  experiments  revealed  a  wealth  of  interesting  nonlinear  phenomena  and  in  particular 
highly  complex  dynamics  including  chaotic  trajectories  as  the  dissipation  parameter  v  decreases. 
Earlier  works  include  the  computations  of  Cohen  et  al.  [4],  Sivashinsky  and  Michelson  [32],  Aimar 
[1],  and  Manneville  [2.3].  Systematic  explorations  of  phase  space  were  carried  out  by  Hyman  and 
Nicolaenko  [15]  and  Hyman  et  al.  [16],  Papageorgiou  and  Smyrlis  [27],  [33]  and  Coward  et  al.  [7] 
who  report  on  many  features  of  the  dynamics.  Kevrekidis  et  al.  [18]  computed  the  bifurcation 
diagram  for  relatively  large  values  of  v,  using  a  bifurcation  algorithm.  For  smaller  values  of  v 
unsteady  phenomena  set  in;  it  was  first  found  by  Papageorgiou  and  Smyrlis  [27]  and  also  Smyrlis 
and  Papageorgiou  [33],  that  in  the  case  of  odd-parity  initial  conditions  there  is  a  period-doubling 
route  to  chaos.  Extensive  numerical  solutions  were  employed  to  provide  strong  evidence  that  the 
route  to  chaos  is  according  to  the  Feigenbaum  scenario  ([8],  [9]);  in  addition  the  universal  constants 
computed  by  Feigenbaum  for  one-dimensional  nonlinear  non-invertible  maps,  were  computed  from 
our  numerical  data  with  two  digit  accuracy.  In  the  sequel  we  term  chaotic  dynamics  just  beyond 
the  accumulation  point,  in  parameter  space,  as  Feigenbaum  chaos. 

In  this  article  we  present  the  results  of  extensive  numerical  computations  for  general  initial 
conditions.  Comparisons  with  previous  studies  (e.g.  [15],  [16],  [18])  have  been  made  where  there  is 
an  overlap  with  full  agreement.  Our  particular  interest  is  in  lower  values  of  v  where  the  dynamics 
gets  increasingly  complicated. 


2  Numerical  methods 

2.1  Approximation  of  the  solution 

The  assumption  that  the  initial  data  is  spatially  periodic  of  period  27r,  allows  us  to  represent  the 
solution  uix.t)  of  the  Kuramoto-Sivashinsky  equation  (from  now  on  KS)  as  : 

uix,t)-  '^(akit)  cos  kx  +  I3kit)  sin  kx)  +  aoii)  (7) 

k>l 
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The  conservative  nature  of  KS  implies  that  ao{t)  remains  constant,  i.e. 

—  — T  /  4“  '^XX  4”  ^'^XXXX^^^  —  0* 

2t\  Jo 

One  could  easily  observe  that  whenever  u{x^t)  is  a  solution  then  so  is  u{x  —  ct^t)  4-  c,  which  allows 
us  to  assume  ao{t)  =  0  for  simplicity.  Replacing  u{x,t)  by  its  Fourier  representation  in  the  PDF 


we  obtain 


Uf  4”  4“  ^xx  4“  ^^xxxx 


where 


=  -  {P  -  vk^)ak  -  Ak)  cos  kx  +  (/?[.  -  {k'^  -  vk‘^)h  -  Bk)  sin  kx) 


^  am/Jn+^  Yl  (oimPn-anPm) 

m-\-n—k  m'~n=k 


Bjz  —  ^  ^  ^  {^m^n  l^m^n)  4“  J  ^  ^  4“  PmPn) 


m-\-n=:k 


m—n=k 


Thus  finally  KS  equation  becomes  equivalent  to  an  infinite  dimensional  system  of  ODEs  ; 

a'k  =  XkOik  +  Ak  (10) 

fii  =  XkPk  +  Bk  (11) 

where  fc  €  IN  and  Xk  =  k"^  -  vk^  with 

Afc  =  Afc(ai,a2,...,/?b/?2,-)  and  Bk  =  Bk{ai,a2,...,^ul32,:.) 

In  the  case  of  the  Generalized  Kuramoto-Sivashinsky  equation 

Uf  “1"  UUx  Uxx  T  ^^xxxx  Bu  0  (1^) 

where  T>  is  a  dispersive  kernel,  see  Papageorgiou  et  al  [25],  the  resulting  infinite  dimensional  system 
of  ODEs  is  quite  similar 

a'f.  =  AfcQfc  +  6l3k  +  Ak  (13) 

^k  -  XkPk  -  ^OLk  +  Bk  (14) 


where  6  is  defined  by 


{Vu){k)  =  iSu{k) 


It  is  already  established  by  H’yashenko  [17],  Collet  et  al  [5],  Goodman  [1-3],  that  all  the  spatial 
Sobolev  norms  of  the  solution  of  KS  equation  for  arbitrary  initial  conditions,  remain  bounded  for 
aU  times,  which  implies  that 

Wfc  =  lim  suplofc^  +  l3k‘^\^^‘^ 


3 


decays  faster  than  any  algebraic  rate.  Collet  et  al.  [6]  have  already  proved  that  the  decay  is 
exponential.  Furthermore  they  present  numericaJ  evidence  to  support  a  conjecture  that  the  number 
of  determining  Fourier  coefficients  is  proportional  to 

Such  numerical  evidence  justifies  approximation  of  the  solution  of  the  KS  by  truncation  of 
higher  frequencies.  The  size  of  the  truncation  can  be  determined  by  the  Collet  et  al.  [6]  numerical 
study  of  the  decay  rate  of  the  Fourier  coefficients,  but  a  first  estimate  of  this  is  obtainable  from 
the  argument  that  follows.  The  KS  equation  in  Fourier  space  can  be  written  as  : 


-u{k,  t)  -  Xku{k,  t)  -  {uu^){k,  t) 


=  Xku{k,t)  -  ik{^u'^){k,t) 


We  also  have  : 


Fact  ;  If  is  positive  and  f{t)  is  in  L°°[0,  +oo)  then  the  solution  of  x'  +  nx  =  f{t)  satisfies: 


lim  sup  |a;(t)|  < -ll/lloo 

++00  fX 


Using  the  above  result  in  combination  with  (17),  we  observe  that,  if  we  had,  even  a  rather  rough 
estimate  of  the  size  of  u  and  its  norms,  we  could  obtain  an  idea  of  the  number  of  Fourier  coefficients 
with  significant  numerical  contribution  to  the  approximate  solution  of  the  KS  equation. 

Having  used  various  numerical  schemes,  we  obtained  sufficiently  good  estimates  of  the  size  of 
u,  as  well  as  of  uux  and  of  \u\iji.  Thus  we  subsequently  obtained  an  estimate  of  the  number  of 
numerically  significant  Fourier  coefficients.  It  should  be  mentioned  here  that  the  existing  analytical 
results  alone  can  not  be  of  much  practical  value  in  determining  the  truncation  size. 

Currently  we  use  a  finite  difference  scheme  of  variable  time  step  to  integrate  the  truncated 
system  of  ODEs 

a'k  =  XkOtk  +  A.k{ct\.,(X2.,  ...,aAr,/3i,/32,  —,fiN)  (19) 

fi'k  =  XkPk  +  Bk{ai,a2,...,aN,fii,P2,—,fiN)  (20) 

where  k  =  1,2, ...,  N  and 

k  k^~^ 

a-mfik-m  +  ~  O-mfim+k)  (21) 


k  k 

Bj^  =  —  ^  ]  i^Ci-ixiOLk—m  fimfik—m)  "I"  'T  ^  '  i,Oim(^m.+k  A  fimfim+k') 


Our  scheme  turns  out  to  be  extremely  stable  due  to  the  form  of  the  hnear  part  of  the  right  hand 
side  of  the  system.  We  emphasize  that  even  though  A*  is  bounded  from  above,  with  a  mild  upper 
bound  namely 
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for  all  A:  €  IN,  it  could  have  a  negative  lower  bound  very  large  in  absolute  in  our  truncated  system. 
For  example,  if  iV  =  32,  i/  =  .01  then  Xn  ^  -9.5  x  10®.  Thus  if  At  is  such  that 

llAflAt  =  0(1) 

then  the  relative  local  numerical  error  for  a^,  /3n  is  also  0(1),  since  at  each  time  step  the  linear 
semigroup  multiplies  oat,  by  exp(/ArAt).  We  should  also  note  here  that  in  the  nonlinear  part  of 
the  system  of  ODEs  corresponding  to  high  frequencies  in  our  truncated  system,  i.e.,  .4^,  most 
of  the  contribution  comes  from  the  low  frequency  coefficients.  This  implies  that  higher  frequency 
coefficients  are  slaved  by  the  low  frequency  ones  due  to  the  scheme;  this  is  a  distinguishing  feature 
of  dissipative  infinite  dimensional  dynamical  systems  and  is  thus  inevitable.  Still  we  wish  to  allow 
the  high  frequencies  as  much  freedom  as  possible  in  order  to  enable  potential  individual  behavior. 
In  achieving  this  we  incorporate  a  modification  suitable  for  stiff  systems  of  ODEs.  Stiffness  is 
caused  here  by  the  big  variation  of  the  values  of  and  thus  the  time  step  in  our  modified  scheme 
varies  according  to  the  size  of  Xk-  Nevertheless  the  relative  local  error  in  every  individual  term  also 
varies,  in  our  scheme,  according  to  its  contribution  to  the  solution. 

Altogether  the  local  numerical  error,  due  to  the  approximation  of  the  solution  of  the  system  of 
ODEs,  is  reduced  to  be  insignificant  compared  to  the  one  caused  by  the  truncation  alone.  It  should 
be  clarified  at  this  point  that  even  with  coarser  time  step  subdivisions  than  the  ones  in  current  use, 
the  relative  local  error  to  the  solution  profile  due  to  the  high  frequencies  is  satisfactory.  High  order 
accuracy,  i.e.  the  local  relative  error  being  kept  below  the  level  of  10~®,  is  employed  only  for  the 
sake  of  studying  high  order  frequencies  and  their  long  time  influence  on  the  nature  of  the  global 
attractor. 


2.2  Numerical  Experiments 


We  have  carried  out  a  very  large  number  of  numerical  experiments,  testing  more  than  800  values 
of  the  dissipation  parameter  i'  ranging  from*  u  =  .99999999  to  v  =  .002.  In  representative  values 
of  our  dissipative  parameter  u  we  carried  out  several  runs  to  make  sure  that  the  behavior  of  the 
attractor  corresponding  to  that  f  does  not  change  with  more  Fourier  coefficients  of  finer  time 
step.  This  is  necessary  since  our  major  concern  is  not  just  the  local  error  but  faithful  reproduction 
of  the  long  time  behavior.  The  size  of  the  truncation  ranged  from  N  —  A  to  N  =  512  making 
computation  in  the  latter  case  rather  slow.  Most  important  of  all  is  that  many  cases  had  to  be 
followed  for  very  long  time  in  order  to  achieve  convergence  to  the  corresponding  attractor  and/or 
to  accurately  classify  the  attractor.  For  example  in  the  case  i/  =  .1212267996068,  the  attractor  is 

*If  I'  >  1  then 


lim  u(x,t)  =uo  =  — 

t— +00  2  7r 


(23) 


since  KS  equation  implies  that 

and  for  u  e  L'^[0,2t],  27r-periodic,  with  u{x)dx  =  0  we  have  the  following  inequalities: 

Wh  <  Wx\2  <  \uxx\2 


On  the  other  hand  if  1/  <  0  the  problem  is  Hi-posed. 


5 


periodic  with  2^^*  distinct  maxima  and  as  many  minima  in  the  norm  as  a  function  of  time.  This 
particular  value  of  u  follows  previous  ones  corresponding  to  13  period  doubling  bifurcations.  The 
period  is  estimated  to  be  T  w  25573,  which  corresponds  to  approximately  5.1  x  10'  time  steps  and 
each  step  corresponding  to  substeps.  Satisfactory  convergence  to  the  periodic  attractor  occurred 
after  10®  time  steps. 

In  representative  cases,  of  different  values  of  v,  we  have  used  our  method  to  estimate  the  error 
caused  by  the  truncation.  The  relative  error,  with  respect  to  the  and  norms  ranged  from 
10“®  to  numerically  insignificant.  Though  10”®  might  look  big,  we  should  point  out  that  our 
study’s  major  target  was  not  to  run  the  most  accurate  algorithm  -  which  we  did  in  selected  cases. 
Instead  it  was  to  carry  out  many  and  long  runs,  and  keep  sufficient  accuracy,  in  order  to  reproduce 
the  nature  of  the  corresponding  attractor.  Classification  of  the  nature  of  the  attractors  at  each 
different  value  of  v  made  it  possible  to  determine  the  windows  in  the  parameter  v  space.  Most  of 
the  runs  had  to  be  carried  out  in  order  to  accurately  determine  the  limits  of  the  windows  of  v,  and 
furthermore  determine  whether  the  transition  between  the  two  kinds  of  attractors  was  smooth  or 
not.  A  case  of  a  smooth  transition  is  a  bifurcation;  A  period-doubling  is  one  such  as  well  as  a  case  of 
eigenvalues  bifurcations  of  the  Jacobian  of  the  flux  of  our  dynamical  system  which  cause  transition 
from  stationary /traveling  attractors  to  periodic  or  chaotic  ones.  A  non-smooth  transition  is  a  case 
when  we  have  for  a  v — interval,  [i^i ,  1^2]  coexistence  of  two  or  more  invariant  locally  attracting 
sets.  While  one  of  them  was  attracting  most  of  the  initial  data,  as  v  varies,  another  one  becomes 
suddenly  more  likely  to  attract  most  initial  data.  For  example,  in  the  case  of  v  -  .232,  we  have 
coexistence  of  a  unimodal  traveling  attractor  and  of  one  with  periodically  appearing  homoclinic 
bursts.  The  second  one  is  more  attracting.  Even  though  the  exact  bifurcation  values  of  p  depend 
on  the  accuracy  of  the  numerical  scheme,  the  nature  of  the  bifurcation  does  not  depend  on  the 
accuracy  provided  that  the  method  is  of  sufficient  accuracy  since  the  system  has  finite  degrees  of 
freedom.  This  observation  forces  us  to  use  higher  accuracy  when  we  seek  better  approximation  of 
the  bifurcation  values  of  p. 


2.3  Characterization  of  the  attractors 


As  p  varies  the  long  time  behavior  of  the  solution  does  also,  in  some  cases  smoothly  and  in  others 
not.  Different  values  of  p  correspond  to  stationary  (unimodal,  bimodal,  trimodal,...),  traveling, 
periodic  and  chaotic  attractors.  Stationary  attractors  are  easy  to  observe,  and  their  stability  can 
be  studied  with  standard  methods.  It  suffices  to  check  the  eigenvalues  of  the  Jacobian  of  the 
truncated  system  at  the  stationary  profile.  Similarly  in  the  case  of  traveling  profiles,  where  one 
could  easily  determine  its  speed  of  propagation  : 


1  0:i/3i  — /?iai  ^  n/A+2^ 


(24) 


where 


=  («i(^  “ 


(CTl!  A)  =  +  ^0) 


In  the  case  of  periodic  attractors  more  sophisticated  numerical  methods  are  required  for  their 
quantification.  First  suspicion  of  a  periodic  attractor  comes  from  the  energy  versus  time  plot  - 
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E{t)  =  \u{t,-)\i2.  The  evidence  becomes  more  convincing  once  we  obtain  the  phase  plane  of  the 
energy  plots  -  i.e.  E{t)  versus  dE{t)fdt,  where  the  values  of  dE{t)fdt  are  accurately  obtainable 
using  a  suitable  interpolation  of  several  points  {tk,E{tk))  -  Periodic  attractors  correspond  to  closed 
phase  curves  and  the  number  of  minima/ma^ma  is  the  number  of  the  points  where  the  phase 
curves  intersect  the  E  -  axis,  i.e.  dE{t)/dt  =  0  with  d?E{t)ldt^  >  0  or  <  0  respectively. 

Nevertheless  we  wish  to  have  more  accurate  quantitative  data  such  as  the  exact  period  as  well  as 
the  exact  values  of  the  local  extrema,  in  order  to  accurately  determine  period  doubling  bifurcations 
and  study  the  potential  fractal  nature  of  the  return  map§  of  the  extrema.  In  order  to  achieve  that, 
we  have  used  a  method  of  calculation  of  the  extrema  by  an  optimized  polynomial  interpolation  over 
a  sufficiently  large  number  of  points  {tk,E{tk))  ■  The  polynomial  of  the  interpolation  is  obtained 
by  a  suitably  weighted  least  squares  method  over  consecutive  points,  many  more  than  the  degree 
of  the  polynomial  sought  for.  Weights  in  the  least  squares  approximation,  reduce  the  effects  of 
higher  order  polynomial  terms  and  the  round-off  error  caused  by  the  computation  of  the  quantity 
used  in  the  interpolation.  It  is  of  critical  importance  that  all  points  of  the  interpolation  lie  in  the 
convex  or  concave  region,  when  we  look  for  minima  or  maxima  respectively.  This  method,  in  the 
case  of  test  functions,  gave  us  accuracy  higher  than  the  machine  precision,  nevertheless,  as  it  could 
be  shown  analytically,  the  precision  in  the  estimation  of  local  extrema,  depend  on  the  time  step 
and  the  degree  of  the  interpolation  polynomial.  Weights  are  introduced  to  make  the  points  which 
are  closer  to  the  local  extremum  have  a  higher  contribution  to  the  interpolation. 

The  original  motive  in  developing  the  interpolation  algorithm  was  to  accurately  determine  the 
period  doubling  bifurcations  and  establish  the  appearance  of  the  Feigenbaum  universal  constants. 
(See  Smyrlis  and  Papageorgiou  [3-3].)  This  method  turned  out  to  be  extremely  helpful  in  deter¬ 
mining  features  such  as  the  fractal  nature  of  attractors,  by  detecting  such  behavior  in  the  return 
map  of  the  energy  minima,  for  example.  Most  successfully  it  has  enabled  us  to  detect  and  classify 
the  quasiperiodicity  of  attractors.  We  should  mention  here  that  the  same  conclusions,  concerning 
periodic  attractors  are  also  obtainable  if  instead  of  the  norm  and  its  extrema  we  had  used  for 
their  detection,  other  quantities  such  as  the  norm  and  the  Fourier  coefficients.  However  there 
is  an  exception  here,  the  case  of  traveling/periodic  attractors,  where  instead  of 

u{x,t  +  T)  =  u{x,t)  for  aU  (a;,t)  €  IR  x  1R+ 


we  have 


u{x,  t  +  T)  =  u{x  -  cT,  t)  for  some  c  €  IR 

Thus  E{t)  has  period  T,  but  u{k,t)  is  in  general  a  quasiperiodic  function  unless  cT /I-k  is  a  rational 
number.  Fourier  coefficients  being  quasiperiodic  produce  very  beautiful  pictures  -  ak(t)  versus 
/?fc(t).  The  speed  of  propagation  and  period  could  also  be  accurately  estimated  using  a  suitable 
interpolation. 

In  the  classification  of  chaotic  attractors  and  transitions  to  chaos  scenarios,  then,  the  extrema 
of  E{t)  and  their  return  maps  played  the  key  role.  We  should  mention  here  that  if  f{t)  is  a 
quasiperiodic  function,  for  example 

f{t)  =  cos  t  +  cos(V^t), 

then,  although  there  is  not  much  to  observe  in  the  graph,  one  could  observe  a  lot  in  the  return 
map.  If  (mjt)fcelN  is  the  sequence  of  the  local  minima,  in  the  order  they  appear  for  t  >  0,  then  their 

^Let  On,  n  €  IN  be  a  sequence  of  real  numbers,  then  the  set  of  points  (on,  On+i),  n  €  IN,  which  is  a  subset  of  IR^, 
is  called  the  Return  Map  of  (onlneiN-  In  the  case  of  a  function  /  :  [0,oo]  ^  IR  the  return  map  of  /  is  the  return  map 
of  the  sequence  of  the  local  minima  (or  maxima  or  extrema)  of  / . 


return  map,  lies  on  a  closed  curve.  In  general  if  f2{t)  are  periodic  functions  with  irrational 

frequency  ratios,  then  the  return  map  of  f{t)  =  fi{t)  +  f2{t)  lies  on  one  or  more  lines.  This  is 
generalized  to  three  or  more  irrational  frequencies. 

Visualization  of  the  return  map  of  the  local  minima  (or  maxima)  of  the  energy  (or  other  quanti¬ 
ties  such  us  the  norm  which  magnifies  the  contribution  of  high  frequencies)  enables  us  to  classify 
the  attractors,  not  only  in  the  quasiperiodic  case,  but  also  in  more  complicated  chaotic  ones. 

In  certain  cases  the  return  map  plot  has  a  fractal  nature.  Foliations  are  observable  when  one 
plots  successive  magnifications  of  the  graphs.  We  are  in  the  process  of  estimating  numerically  the 
HausdorfF  dimension  of  the  graph  corresponding  to  the  return  map  in  order  to  establish  the  rate  of 
growth  of  the  HausdorflF  dimension  as  i/  decreases,  particularly  in  those  chaotic  cases  where  there  is 
no  recognizable  pattern.  In  such  studies  it  is  imperative  that  the  local  extrema  are  computed  with 
a  high  accuracy,  in  particular  if  magnifications  in  the  return  map  of  several  orders  of  magnitude 
are  to  be  carried  out. 

3  Numerical  results 

3.1  Dynamics  for  “large”  i/t  0.09  <  u  <1 

For  values  of  i/  larger  than  about  0.1  a  fairly  complete  picture  of  the  dynamics  has  been  given 
by  numerous  authors  ([15],  [16],  [18]).  Steady  state  or  traveling  waves,  including  their  bifurcation 
diagrams  are  computable  by  well-developed  bifurcation  algorithms,  see  for  example  ([18]).  Our 
interest  is  on  the  spatio-temporal  evolution  in  chaotic  attractors  which  led  to  the  development  of 
the  accurate  code  tracking  the  time  evolution  of  solutions,  described  in  Section  2.  For  completeness, 
however,  and  to  set  the  stage  for  the  results  that  follow  we  give  a  brief  description  of  the  dynamics 
corresponding  to  “large”  i/.  We  note  that  as  t/  decreases  computed  attractors  are  not  necessarily 
global  ones;  our  numerical  results  are  based  on  large  time  solutions  of  general  initial  value  problems 
and  the  numerical  solutions  given  here  belong  to  the  most  strongly  attracting  ones  in  the  sense 
that  they  are  reproduced  by  fairly  different  initial  conditions  (e.g.  a  sinusoidal  initial  condition  or 
one  having  a  white  noise  spectrum). 

A  schematic  of  the  dominant  attractors  in  this  range  is  given  in  Figure  1.  The  interval  v  >1 
is  not  included  in  Figure  1  since  it  is  trivial  in  the  sense  that  a  uniform  steady  state  is  achieved 
here  which  is  equal  to  the  spatial  mean  of  the  initial  data.  The  figure  is  drawn  to  scale  in  order 
to  give  the  relative  sizes  of  attractors.  The  letters  A-F  are  used  to  identify  the  different  inter¬ 
vals  which  are  briefly  described  below.  According  to  extensive  computations,  there  is  an  interval 
(.12116,  .1212267996064)  which  cannot  be  distinguished  on  Figure  1  due  to  its  small  size  -  we  call 
this  interval  Ei  and  describe  it  below. 

Interval  A:  0.307602  <v<\. 

For  0.30765  <  v  <  \  solutions  get  attracted  to  a  27r-periodic  steady  state  (unimodal  fixed  point). 
The  energy,  or  the  spatial  X^-norm  of  the  solution,  defined  as  E{t)  —  v?{x,t)dx  increases  as  a 
function  of  u  in  this  interval  and  attains  a  maximum  at  //  «  0.59094. 

Interval  B:  0.232491  <v<  0.307601. 

For  0.232491  <  v  <  .307601  a  27r-periodic  (unimodal)  steady-state  traveling  wave  is  found  -  the 
wavespeeds  increase  monotonicaUy  from  zero  near  the  upper  edge  of  the  window  to  approximately 
0.956  at  1/  =  0.232491. 
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Interval  C:  0.17735  <v<  0.23249. 

A  bifurcation  occurs  just  below  v  =  0.232491,  and  in  the  range  0.17735  <v<  0.23249  periodic 
homoclinic  bursts  (PHB)  emerge  as  the  most  attracting  states.  The  PHB  solutions  found  in  this 
window  have  X^-norms  which  are  constant  except  for  periodically  occurring  homoclinic  bursts.  The 
time  intervals  between  bursts  vary  from  approximately  21  -  30  time  units  for  v  between  0.23249 
and  0.2,  and  then  increase  dramatically  as  v  decreases  further;  for  example  the  time  between  bursts 
for  V  =  0.179,0.178,0.1775,0.1774,0.177.35  are  approximately  200,450,1700,4500  and  18000  time 
units  respectively. 

Representative  results  are  given  in  Figures  2  and  3  for  =  0.2.  Figure  2  shows  the  evolution 
of  E{t)  after  the  solution  has  settled  into  an  attractor  -  a  range  between  750  and  1000  time  units 
is  shown.  The  energy  is  constant  with  homoclinic  bursts  occurring  at  approximately  every  30  time 
units.  The  phase  plane  of  E{t)  provides  strong  evidence  that  the  bursts  are  periodic.  The  solution 
spends  most  of  its  time  at  the  rightmost  point  in  the  phase  plane  diagram  and  every  30  time  units  it 
loops  around  in  less  than  5  time  units,  with  the  process  repeating  periodically.  The  profile  between 
bursts  is  steady  for  30  time  units  and  is  in  fact  a  bimodal  fixed  point,  that  is  only  even  Fourier 
modes  are  non-zero.  It  has  been  confirmed  that  the  bimodal  profile  during  bursts  is  obtainable 
from  the  unimodal  global  fixed  point  attractor  at  i/  =  0.8  by  use  of  the  property  (26)  discussed 
in  Section  2.  We  can  surmise,  therefore,  that  the  bimodal  fixed  point  attractor  is  unstable  and 
overlaps  with  a  time-periodic  attractor. 

The  evolution  of  E{t)  does  not  provide  details  about  individual  Fourier  coefficients  since  it  is  an 
integral  quantity.  The  time  history  of  the  first  two  harmonics  is  given  in  Figure  3.  The  coefficients 
of  cos(a;),  sin(x)  and  cos(2a;),  sin(2a;)  are  plotted  on  the  same  scale  so  that  their  evolution  can  be 
followed  concurrently.  It  is  evident  from  the  primary  Fourier  mode  plots  that  the  amplitudes  are 
zero  except  during  bursts  which  appear  as  sharp  peaks.  The  peaks  alternate  from  small  size  to  larger 
size  during  each  burst;  there  also  seems  to  be  an  energy  exchange  between  the  two  amplitudes  in 
that  during  each  burst  the  cos(a;)  and  sin(a;)  amplitudes  attain  alternatively  large  and  small  peaks. 
The  sign  of  individual  amplitudes  during  bursts  does  not  appear  to  have  a  recognizable  pattern 
but  seems  to  be  random. 

The  situation  is  much  more  regular  for  the  2nd  Fourier  mode.  The  time  evolution  of  the 
amplitudes  of  cos(2a;)  and  sin(2a;)  are  now  in  phase  with  bursts  appearing  as  sharp  jumps  connecting 
equal  but  opposite  steady  states.  Similar  results  persist  for  higher  modes,  the  main  difference  being 
a  marked  decrease  in  amplitude. 

Interval  D:  0.131815  <u<  0.1773. 

The  time  increase  between  bursts  heralds  the  appearance  of  a  new  fixed  point  attractor.  It  is 
found  that  in  the  range  0.131815  <  v  <  0.1773  solutions  get  attracted  to  bimodal  steady  states, 
i.e.  TT-periodic  stationary  profiles.  According  to  the  property  (26)  these  solutions  are  obtainable 
from  the  unimodal  solutions  in  the  subinterval  [0.5272, 0.7092]  of  A  by  choosing  p  =  2. 

Intervals  E,  Ei:  0.1212267996068  <u<  0.13181,  and  £3:  0.12116  <u<  0.1212267996064. 

The  bimodal  steady  states  loose  stability  and  a  time  periodic  attractor  is  entered.  The  energy 
E(t)  has  an  approximate  period  of  0.93862  time  units  at  the  beginning  of  the  window,  v  =  0.13181. 
The  period  increases  monotonically  as  u  is  decreased  and  the  first  period  doubling  takes  place  in 
the  interval  .12175  <  v  <  .121752.  This  pattern  is  repeated,  i.e.  a  monotonic  increase  in  the 
period  and  then  a  period  doubling,  and  the  following  eleven  successive  period  doublings  were  found 
in  each  of  the  open  intervals  (.1213145,  .121315),  (.12124493,  .12124495),  (.1212.3065,  .1212307), 
(.12122762,  .12122775),  (.121226975,  .12122698),  (.121226835,  .12122684),  (.121226805,  .121226808), 
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(.1212268012,  .1212268014),  (.1212267998,  .1212268),  (.12122679968,  .1212267997)  and 
(.121226799622,. 121226799625).  The  period  by  the  twelveth  period  doubling  is  T  «  12780  time 
units.  Figure  4  shows  the  phase  planes  of  the  energy  for  representative  cases  capturing  the  first 
five  period-doublings.  The  values  of  v  are  given  on  the  Figure  as  well  as  the  number  of  minima 
of  the  energy  -  there  are  2^  minima  where  p  is  the  number  of  period  doublings  that  have  taken 
place.  It  has  been  confirmed,  but  is  not  included  here,  that  the  cascade  follows  the  Feigenbaum 
scenario  ([8],  [9])  with  r'-intervals  between  successive  period  doublings  decreasing  geometrically  by 
a  universal  factor  6  =  4.6692016 . . ..  The  data  analysis  and  methods  used  are  identical  to  those 
of  Papageorgiou  and  Smyrlis  [27]  and  Smyrlis  and  Papageorgiou  [33]  who  considered  odd-parity 
solutions  of  the  KS  equation.  The  second  Feigenbaum  universal  constant  which  pertains  to  the 
geometry  of  the  attractor  is  also  accurately  supported  by  those  results.  By  analogy  with  the  theory 
of  non-invertible  one- dimensional  maps,  coupled  with  the  numerical  evidence  given  here,  it  is  rea¬ 
sonable  to  assume  that  an  accumulation  point  exists,  vq  say,  below  which  the  solution  is  chaotic; 
this  of  course  is  hard  to  prove  since  we  are  dealing  with  differential  systems  of  large  dimensions  (24 
coupled  equations  are  used  in  this  window). 

After  extensive  computations  we  are  able  to  identify  an  interval  Ei  =  [0.12116, 0.1212267996064] 
which  we  describe  as  supporting  Feigenbaum  chaos.  Evidence  for  this  comes  from  the  study  of 
return  maps  of  the  energy  minima  (see  Section  2  for  a  description  of  the  data  analysis  tools  used 
here)  which  give  “strange  attractors”  with  as  many  as  2^^  pieces  at  the  beginning  of  interval  E\. 
The  number  of  pieces  decrease  and  when  v  =  0.12118  the  attractor  consists  of  only  one  piece.  A 
picture  of  a  strange  attractor  constructed  from  the  return  map,  is  given  in  Figure  5(i)-(vii)  which 
corresponds  to  v  =  0.12122679907.  The  return  map  of  a  large  number  of  energy  minima  beyond 
transient  behavior  is  shown  in  Figure  5(i).  There  are  32,034  points  in  this  plot  and  each  point  is 
represented  by  a  dot.  The  points  lie  on  a  regular  looking  object  and  since  each  point  gets  mapped 
to  another  point  in  the  plot  by  the  action  of  the  dynamical  system,  the  pieces  traced  out  mimic  the 
analogous  one-dimensional  non-invertible  map  representation  of  the  flow;  the  equivalent  plot  for 
the  logistic  map  would  trace  out  parts  of  a  quadratic  function.  Interestingly,  Figure  5(i)  resembles 
the  quadratic  map  but  the  differences  have  not  been  quantified.  The  similarity  property  of  strange 
attractors  is  however  exhibited  in  Figure  5.  Starting  from  Figure  5(i),  successive  enlargements  are 
made  of  the  pieces  of  the  attractor  which  are  marked  by  the  letter  ‘e’.  The  scale  of  the  enlargement 
can  be  followed  by  noting  the  horizontal  scale  of  successive  diagrams.  The  first  enlargement.  Figure 
5(ii),  is  of  what  appear  as  two  pieces  in  5(i)  near  the  minimum  of  the  plot.  On  enlargement  a  picture 
which  is  clearly  similar  to  the  complete  attractor  in  5(i)  emerges.  The  two  pieces  near  ‘e’  in  5(ii) 
are  enlarged  to  give  5(iii);  again  similarity  with  the  previous  plots  is  remarkable  considering  that  an 
amplification  of  approximately  a  factor  of  50  has  taken  place.  The  one  piece  near  ‘e’  in  5(iii)  is  now 
enlarged  to  give  the  three  clusters  in  5(iv).  The  rightmost  marked  piece  is  enlarged  to  give  5(v)  and 
its  rightmost  piece  is  enlarged  to  give  5(vi).  Finally  an  enlargement  of  the  leftmost  marked  cluster 
of  5(vi)  gives  the  final  plot  Figure  5(vii)  which  contains  just  four  points.  The  successive  plots  5(iv)- 
(vii)  provide  additional  strong  evidence  of  the  self- similarity  of  the  attractor.  We  note  that  such 
conclusions  are  only  possible  if  a  large  and  accurate  data  base  is  available  -  for  instance  in  going 
from  Figure  5(i)  to  Figure  5(vii),  an  amplification  factor  of  approximately  3  x  10®  has  taken  place 
along  the  horizontal  axis.  In  the  absence  of  a  large  and  accurate  data  set  of  minima  (see  Section  2 
for  the  interpolation  algorithms),  the  number  of  return  map  “iterations”  which  produce  the  strange 
attractor  would  be  inadequate  for  sufficient  levels  of  enlargement  that  exhibit  self-similarity.  In 
addition  loss  of  accuracy  in  the  minima  estimation  would  produce  a  noisy  data  set  which  is  of  little 
value  in  exhibiting  self-similarity. 
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Interval  F:  0.09051  <  <  0.121155. 


Near  the  upper  edge  of  this  window,  0.121155  <  i'  <  0.12105,  solutions  are  attracted  to  ho¬ 
moclinic  bursts,  with  the  bursts  being  chaotic  but  with  a  low  dimensional  attractor.  The  bursts 
are  spaced  apart  at  roughly  equal  times  and  in  between  them  the  profile  is  bimodal  and  is  obtain¬ 
able  from  the  equivalent  unimodal  steady  state  in  interval  A  according  to  the  self-similar  property 
(26).  Typical  results  are  given  for  v  =  0.12115.  The  homochnic  bursting  is  shown  in  the  evo¬ 
lution  of  the  energy  and  its  corresponding  phase-plane.  Figures  6;  the  phase-plane  indicates  the 
non-periodicity  of  the  bursts.  It  has  been  confirmed  that  the  profile  in  between  bursts  is  given  by 
u(a;;  0.12115)  =  217(2x;  0.4846)  where  U{y;iy),  0  <  y  <  27r  is  the  steady-state  profile  in  interval  A. 
Further  information  about  the  chaotic  bursts  is  given  in  Figures  7(a)-(b).  Figure  7(a)  gives  the 
evolution  with  time  of  the  critical  points  of  E{t),  that  is  its  local  maxima  and  local  minima.  For 
a  given  time  interval,  there  is  a  discrete  set  of  energy  maxima  and  minima  and  each  of  these  is 
represented  by  a  dot  in  the  plot.  The  longer  and  relatively  regularly  spaced  gaps  with  no  points 
in  the  plot  denote  time  intervals  between  bursts  while  darker  segments  correspond  to  the  bursts. 
The  irregularity  of  the  positioning  of  the  points  indicates  the  chaotic  aspect  of  the  bursts  and  this 
can  be  seen  more  clearly  in  the  return  map  of  the  energy  minima  in  Figure  7(b)  which  suggests 
a  strange  attractor  with  foliations.  The  bimodal  steady  state  in-between  bursts  is  at  the  point 
(10.7415876, 10.7415876)  on  the  return  map  diagram,  which  is  marked  with  an  open  circle.  Succes¬ 
sive  enlargements  of  the  attractor  are  given  in  Figure  8(i)-(vi).  Again  the  point  corresponding  to 
the  inter-bursting  bimodal  steady  state  is  marked  with  an  open  circle  and  five  successive  enlarge¬ 
ments  are  made  in  its  neighborhood.  The  similarity  of  the  attractor  is  discernible,  particularly  if 
we  note  the  horizontal  extent  of  the  attractor  given  in  each  plot  in  Figures  (ii)-(vi);  the  horizontal 
extents  are  approximately  0.6,  0.15,  0.011,  0.001  and  0.000005  for  Figures  8(ii)-(vi)  respectively. 
An  eidargement  of  the  order  10^  has  taken  place  in  going  from  Figure  8(ii)  to  Figure  8(vi)  and 
the  positioning  of  the  dots,  that  is  the  geometry  of  the  attractor,  is  remarkably  similar;  from  a 
numerical  point  of  view,  we  note  that  such  details  can  only  be  captured  if  the  energy  minima  are 
computed  accurately  to  at  least  seven  decimal  places  -  this  is  another  reason  for  the  high  order 
interpolation  schemes  used  in  constructing  the  minima  from  a  set  of  energy  values  at  discrete  time 
intervals  no  smaller  than  the  time-step  of  the  integrator. 

As  the  value  of  v  is  decreased  below  approximately  0.121045  the  homoclinic  bursts  take  place 
more  regularly  in  time  and  they  are  almost  periodic.  A  more  careful  scrutiny  of  the  data  from  the 
representative  case  v  =  0.1  reveals  that  the  flow  is  chaotic,  and  the  bursting  phenomenon  has  a  low 
dimensional  attractor.  Once  again  the  unstable  quiescent  states  between  bursts  give  profiles  which 
are  bimodal  and  given  by  u(a;;0.1)  =  2U{2x;0.4).  A  plot  of  E{t)  together  with  a  magnification  of 
a  burst  is  shown  in  Figure  9.  The  bursts  are  almost  periodic  and  take  place  every  33  time  units, 
approximately;  the  energy  phase-plane  has  been  shown  not  to  be  a  closed  curve  of  index  2  as  in 
the  homochnic  burst  solutions  found  in  interval  C  described  above.  Figure  10  shows  the  return 
map  obtained  from  the  energy  minima.  The  point  marked  P  corresponds  to  the  steady  state  in 
between  bursts  and  as  the  magnification  shows,  the  attractor  in  the  neighborhood  of  this  fixed 
point  consists  mainly  of  two  fines  meeting  at  an  angle  smaller  than  7r/2.  Additional  information 
about  the  bursting  is  given  in  Figure  11.  The  Fourier  amplitudes  of  the  first  four  inodes  are  plotted 
against  time,  that  is  we  plot  -f  6?  for  i  =  1,4  where  u{x,t)  =  E^_i(a„  cos(nar)  -I-  sin(na;)) 
is  the  Fourier-Galerkin  representation  of  the  solution.  In  between  bursts,  only  the  even  modes  are 
non-zero  and  as  seen  from  the  plots  the  amplitude  of  the  fourth  mode  is  smaller  than  that  of  the 
second  one.  During  a  burst,  however,  the  even  modes  loose  energy  to  the  odd  modes  as  can  be 
seen  clearly  from  the  Figure;  the  amplitude  of  even  modes  always  decreases  (except  right  at  the 
end  of  the  burst)  and  at  the  same  time  the  amplitudes  of  odd  modes  become  non-zero  for  the 
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duration  of  the  burst  but  vanish  again  after  the  homoclinic  burst  is  finished.  This  behavior  is  very 
characteristic  of  the  global  energy  transfer  mechanism  present  in  the  KS  equation. 

3.2  Dynamics  for  lower  values  of  v  <  0.09. 

As  V  decreases  the  dynamics  become  more  complicated.  A  larger  number  of  determining  modes  is 
required  and  co-existence  of  attractors  is  more  likely.  From  a  numerical  viewpoint  the  computations 
become  more  expensive  and  accuracy  criteria  more  stringent. 

A  large  number  of  numerical  experiments  have  been  carried  out  over  several  years  and  a  sum¬ 
mary  of  the  attractors  is  given  in  Figure  12.  The  Figure  represents  the  z/-line  (not  drawn  to  scale) 
along  with  the  most  strongly  attracting  large  time  solutions  found  by  us;  note  that  the  various 
attractors  given  in  Figure  12  are  not  necessarily  global  ones  but  the  window  boundaries  given 
delineate  the  regions  where  an  attractor  loses  stability.  The  type  of  attractor  was  determined  by 
applying  the  data  analysis  described  earlier. 

In  interval  E  described  earlier,  we  provided  evidence  of  a  period-doubling  route  to  chaos,  and 
also  a  route  to  chaos  through  homoclinic  bursting  as  in  interval  F  above.  A  third  route  to  chaos 
is  via  a  quasi-periodic  attractor  and  this  route  is  supported  by  the  first  three  windows  in  Figure 
12.  A  time  periodic  attractor  looses  stability  to  a  quasi-periodic  one  dA  v  =  0.087679  and  there 
is  transition  to  chaos  below  v  -  0.08702.3.  It  is  quite  hard  to  discern  the  difference  between  a 
chaotic  attractor  and  a  quasi-periodic  one  by  inspection  of  the  signal  of  E{t).  The  return  map  is 
a  valuable  tool  in  doing  this,  since  the  return  map  (of  the  energy  minima  for  instance)  of  a  quasi 
periodic  attractor  consists  of  lines  with  no  foliations;  an  analogy  is  the  return  map  obtained  from  a 
Poincare  section  of  motion  on  a  torus,  in  which  case  it  is  a  circle  whose  perimeter  becomes  dense  as 
the  iterations  increase.  In  the  chaotic  regime,  the  return  map  contains  foldings  and  self-similarity 
as  exhibited  in  several  examples  earlier.  This  type  of  chaotic  motion  is  differentiated  from  vrhat 
we  term  as  unrecognizable  chaos  in  that  the  tools  employed  by  us  to  analyze  the  dynamics,  yield 
information  about  the  structure  of  the  attractor. 

As  can  be  seen  from  Figure  12,  among  the  regions  of  different  unsteady  and  complicated  spatio- 
temporal  motions  there  emerge  stable  fixed  point  attractors  with  increasingly  higher  modal  be¬ 
havior,  i.e.  with  shorter  periods.  In  Figure  13  we  give  a  summary  of  the  computed  attractors  so 
far,  as  v  decreases  to  approximately  0.006.  The  numbers  on  the  figure  denote  the  modal  form  of 
the  solution,  for  example  3M  is  a  trimodal  attractor  with  the  only  non-zero  Fourier  coefficients 
being  multiples  of  three.  As  already  noted  in  the  Appendix,  one  can  construct  steady  states  of 
increasingly  shorter  periods  (i.e.  higher  modal  behavior)  by  application  of  the  similarity  property 
(26).  This  construction  does  not  imply  stability  but  our  numerical  results  produced  by  solving 
initial  value  problems,  show  that  such  fixed  points  are  observable  as  large  time  solutions,  given 
appropriate  initial  conditions.  In  what  follows  we  use  numerically  computed  solutions  from  each 
of  the  windows  which  support  multimodal  steady  states  (i.e.  3M,  4M,  5M,  6M,  7M,  8M  and  9M) 
and  show  that  they  come  from  the  unimodal  fixed  point  attractor  in  interval  A  (see  Figure  1). 

Making  the  assumption  that  the  multimodal  fixed  points  arise  from  the  self-similar  property 
(26),  we  first  transform  the  i/-windows  which  support  these  steady  states  by  multiplying  by  iV^ 
where  N  is  the  modal  characterization  of  the  attractor.  This  transformation  gives  the  corresponding 
viscosity  interval  in  the  global  fixed  point  attractor  A  of  Figure  1.  The  results  are  summarized  in 
Table  1. 
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Attractor 

1/  interval 

Transformed  zz-interval  (xiV^^) 

3M 

[.07465, .059922] 

[.67185,.539298] 

4M 

[.043848,.034625] 

[.701568,.554] 

5M 

[.026984,.022604] 

[.6746,.5651] 

6M 

[.01829,.015845] 

[.65844,.57042] 

7M 

[.012994,.01174] 

[.636706,.57.526] 

8M 

[.009905,.009204] 

[.63.392,.589056] 

9M 

[.007479, .007.305] 

[.605799,.59170.5] 

lOM 

.006 

.6 

Table  1 

Table  1  indicates  that  the  transformed  intervals  for  each  of  the  multimodal  attractors  listed  lie  in 
the  i/-interval  which  supports  a  unimodal  global  fixed  point.  In  what  follows  we  provide  numerical 
evidence  that  each  of  the  multimodal  attractors  at  decreasing  values  of  v,  are  indeed  obtainable 
from  the  transformation  (26).  This  is  done  for  the  profile  corresponding  to  v  —  .6  in  interval  A,  since 
this  is  seen  to  be  common  to  all  entries  in  the  Table.  For  iV  =  3, . . . ,  10  we  compute  the  steady  state 
attractor  at  the  corresponding  values  v  —  0.6/iV^  and  denote  the  profile  by  un{x).  These  profiles 
are  shown  in  Figure  14;  the  period  of  corresponding  profiles  is  seen  to  be  2Tr/N  and  the  amplitude 
increases  with  N.  Note  that  due  to  translation  invariance  each  of  the  profiles  can  be  shifted  by  a 
constant  along  the  a;-axis  and  it  remains  a  periodic  solution.  This  observation  is  particularly  useful 
in  the  normalization  construction  that  follows.  Each  profile  is  then  normalized  according  to  (26), 
that  is  we  define  a  new  function,  17jv(?/)  say,  on  [0,27r]  by  applying  the  transformation 

UN{y)  =  (25) 

This  transformation  was  first  noted  by  Frisch  et  al.  [12]  who  went  on  to  carry  out  the  linear 
stability  of  these  profiles  for  large  N  -  see  below.  If  the  self-similarity  property  (26)  is  valid  for 
these  multimodal  fixed  points,  then  the  UN{y)  for  iV  =  1, . . .,  10  are  identical  to  within  a  horizontal 
shift  and  equal  to  the  profile  obtained  at  i/  =  0.6.  It  was  found  that  the  computed  profiles  4M-10M 
had  to  be  shifted  by  w  in  order  to  put  them  in  phase  with  the  globally  attracting  profile  at  =  0.6 
(see  earlier  comments).  The  results  are  given  in  Figure  15  with  different  symbols  used  to  denote 
different  values  of  N.  The  solid  line  corresponds  to  the  generating  fixed  point  at  z/  =  0.6.  It  is 
seen  from  the  Figure  that  the  self-similarity  property  is  indeed  the  controlling  mechanism  in  the 
generation  of  these  high-modal  steady  states  as  u  decreases. 

The  data  of  Table  1  indicate  the  parameter  range  of  the  fixed  point  unimodal  attractor  which 
needs  to  be  considered  in  order  to  transform,  according  to  the  similarity  property,  to  the  interval 
obtained  numerically  for  a  given  multimodal  attractor.  In  the  stability  studies  of  Frisch  et  al  [12] 
and  Papageorgiou  et  al  [26],  the  underlying  basic  profiles  are  the  multimodal  ones  indicated  above 
and  constructed  in  the  Appendix,  as  the  value  of  iV  or  p  increases.  Such  solutions  are  found  to  be 
stable  in  the  z/-interval  (0.6, 0.7),  approximately.  This  result  would  imply  that  given  a  value  of  u 
in  the  stability  window,  the  multimodal  profiles  obtained  by  the  self  similarity  property  (26),  i.e. 
the  one-parameter  family  Nu{Nx;  i>lN^)  for  iV  =  1, 2, . . .  as  iV  gets  large,  are  linearly  stable.  Our 
numerical  results  are  in  agreement  with  this  theoretical  prediction  for  N  as  large  as  10.  Note  that 
computation  of  steady  states  with  iV  >  10  has  not  been  achieved  yet  because  the  window  which 
supports  such  states  has  size  of  the  order  of  machine  precision.  This  may  be  achievable  with  higher 
precision  arithmetic  but  such  explorations  are  not  undertaken  here. 
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4  Conclusions 


We  have  given  many  quantitative  features  of  solutions  of  the  Kuramoto-Sivashinsky  equations 
computed  with  general  initial  data  and  spatially  periodic  boundary  conditions.  Features  such  as 
strange  attractors  and  periodic  or  chaotic  bursting  phenomena  have  been  elucidated.  Of  particular 
interest  is  a  set  of  computed  multimodal  solutions  found  at  decreasing  values  of  u  and  supported 
on  smaller  and  smaller  windows.  It  has  been  shown  that  these  solutions  (the  last  computed  one  is 
a  decamodal  profile,  i.e.  one  with  aU  Fourier  coetficients  zero  except  those  which  are  multiples  of 
10)  derive  from  the  unimodal  fixed  point  attractor  in  .3  <  <  1,  by  a  self-similar  transformation 

property  of  the  equation.  It  can  be  concluded  that  such  solutions  are  at  least  linearly  stable  (in 
the  sense  that  they  are  computable  as  large  time  solutions  of  a  nonlinear  initial  value  problem) 
and  our  numerical  results  are  in  full  agreement  with  the  stability  theory  of  Frisch  et  al  [12]  and 
the  modulation  theory  of  Papageorgiou  et  al  [26]  who  study  the  stability  of  such  multimodal 
steady  states.  The  numerical  results  given  here  provide  strong  support  for  the  modulation  stability 
theories. 


APPENDIX 

The  following  property  of  (2)  is  useful.  It  has  been  established  numerically  that  for  0.307602  < 
V  <l  (this  is  interval  A  described  in  Figure  1  and  Section  3)  the  KS  has  a  global  unimodal  fixed 
point  attractor.  That  is  given  fairly  general  initial  conditions  (for  example  a  white  noise  spectrum) 
solutions  evolve  to  a  steady  state  with  period  2n.  Denote  such  steady  states  by  U{x\  u)  with  v  in 
the  interval  A.  For  any  constants  p  and  c,  equation  (2)  has  the  following  two-parameter  family  of 
solutions 

u{x,t)  =  pU{p(x  -  cty, p^v)  +  c.  (26) 

This  follows  by  the  scale  and  Gahlean  invariance  of  KS.  For  stationary  waves  c  =  0,  while  peri¬ 
odicity  implies  that  p  is  an  integer.  Property  (26)  can  be  used  to  generate  steady  states  of  (2) 
at  geometrically  decreasing  values  of  the  viscosity.  For  example  if  the  27r-periodic  steady  state  is 
known  at  i/  =  0.8,  we  can  construct  a  steady  state  at  z/  =  0.8/2^  by  choosing  /)  =  2;  this  solution 
has  period  tt,  i.e.  is  bimodal,  and  ampUtude  twice  that  of  the  starting  steady  state;  this  folding  and 
scaling  process  can  be  repeated  ad  infinitum  to  obtain  multimodal  steady  states  at  geometrically 
decreasing  values  of  v.  Clearly  in  the  small  v  limit  the  amplitude  of  such  steady  states  is  of  the 
order  of  which  is  also  the  order  of  the  X^-norm  of  these  steady  state  solutions.  Even  though 

the  steady  states  exist  and  can  be  constructed  as  shown  above,  they  are  not  necessarily  stable  and 
may  not  be  observed  in  numerical  experiments.  What  has  been  found,  however,  is  that  there  is  a 
subwindow  of  interval  A  which  can  be  used  to  generate  multimodal  steady  states  which  are  linearly 
stable. 
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Figure  1  Schematic  of  the  various  attractors  for  larger  values  of  v.  The  diagram  is  to  scale 
and  the  various  most  strongly  attracting  solution  regions  are  indicated. 


17 


sin  cos  sin  cos 


750  800  850  900  950  1000 

t 


Figure  3  Interval  C:  Periodic  homoclinic  bursts  for  i/  =  0.2.  Time  evolution  of  the  1st 
and  2nd  Fourier  modes.  The  graphs  show  the  evolution  of  the  coefficients  of  cos(a;),  sm(a;)  and 
cos(2a;),sin(2x).  Higher  mode  amplitudes  have  much  lower  amplitudes. 
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Figure  4  Interval  E:  Period  doubling  route  to  chaos  according  to  a  Feigenbaum  scenario.  The 
Figure  gives  the  phase  planes  of  the  energy  as  v  decreases.  The  values  of  v  and  the  number  of 
relative  minima  in  each  curve  are  indicated  on  the  Figure.  Chaos  sets  in  as  the  phase  plane  attains 
increasingly  more  turns  which  become  infinitely  many  as  an  accumulation  point  is  reached. 
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Figure  5  Interval  Ei:  Feigenbaum  chaos,  u  -  0.12118.  Strange  attractor  produced  by  the 
return  map  of  the  energy  minima.  The  Figure  gives  successive  enlargements  near  the  point  ‘e’ 
marked  on  each  subplot  indicating  the  self-similar  nature  of  the  attractor.  An  overall  enlargement 
of  3  X  10®  takes  place  in  going  from  5(i)  to  5(vii). 
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Figure  6  Interval  F:  Chaotic  homoclinic  bursts,  v  =  0.12115.  The  energy  and  its  phase  plane 
are  shown.  The  bursts  occur  at  regular  intervals  and  are  chaotic  in  nature.  The  steady  attractor 
between  bursts  corresponds  to  the  unimodal  global  steady  attractor  in  interval  A  having  u  =  0.4846 
and  the  similarity  transformation  (26)  (see  text). 
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Figure  7  Interval  F:  Chaotic  homoclinic  bursts,  u  =  0.12115.  Time  history  of  the  critical  points 
of  the  energy  (both  maxima  and  minima).  The  absence  of  points  indicates  a  constant  energy.  An 
extended  chaotic  region  is  seen  between  approximately  t  =  1000  -  1400,  which  then  gets  attracted 
to  chaotic  homoclinic  bursts.  The  lower  figure  shows  the  return  map  of  the  energy  minima.  A 
strange  attractor  emerges  and  the  point  between  bursts  is  indicated  by  an  open  circle. 
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Figure  8  Interval  F:  Chaotic  homoclinic  bursts,  u  =  0.12115.  Self-similarity  of  the  strange 
attractor.  Successive  enlargements  are  shown  in  the  neighborhood  of  ‘o’.  An  overall  enlargement 
of  the  order  of  10'^  takes  place  in  going  from  8(i)  to  8(vi). 
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Figure  10  Interval  F:  Chaotic  homoclinic  bursts,  u  =  0.1.  Return  map  of  the  energy  minima 
indicating  a  strange  attractor.  The  point  P  indicates  the  position  in  between  bursts. 
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Figure  11  Interval  F:  Chaotic  homoclinic  bursts,  u  =  0.1.  Evolution  of  the  first  4  Fourier 
amplitudes  together  with  the  energy  (sum  of  the  Fourier  amplitudes).  It  is  seen  that  in  between 
bursts  the  even  Fourier  modes  are  non-zero  while  the  odd  ones  are  zero,  and  the  steady  attractor 
is  bimodal,  then.  During  bursting  an  energy  exchange  takes  place  with  odd  modes  gaining  energy 
from  even  modes.  The  interval  between  bursts  remains  fairly  constant  over  long  time  periods. 
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Figure  12  Schematic  of  the  various  computed  attractors  for  smaller  values  of  nu  {v  <  0.09). 
The  diagram  is  not  drawn  to  scale  but  window  boundaries  are  given  along  with  a  brief  description 
of  the  type  of  computed  dynamics.  The  attractors  are  not  global. 
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Figure  13  Distribution  of  computed  steady  attractors  for  v  <  0.09.  The  diagram,  drawn  to 
scale,  shows  the  z/-intervals  where  steady  state  solutions  were  found.  The  notation  IM,  3M,  etc. 
denotes  profiles  which  axe  unimodal,  trimodal  etc..  The  windows  of  support  of  higher  modal  steady 
states  get  decreasingly  small  as  v  decreases  and  their  computation  becomes  prohibitively  difficult. 
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Figure  14  Computed  profiles  corresponding  to  values  of  v  taken  from  the  global  attractor 
for  large  v  (interval  A  in  Figure  1)  and  the  intervals  3M,  4M,  5M,  6M,  7M,  8M,  9M  and  lOM. 
The  corresponding  values  of  nu  are  indicated  on  the  Figures  and  are  chosen  in  order  to  check  the 
self- similarity  of  the  solutions  (see  text). 
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Figure  15  Similarity  of  steady  states  as  v  decreases.  Collapse  of  the  profiles  given  in  Figure 
14  onto  the  universal  profile  corresponding  to  v  =  0.6.  Different  symbols  correspond  to  data  taken 
from  the  indicated  intervals.  The  similarity  property  (26)  is  used  in  the  construction  of  the  Figure 
-  see  text  for  details. 
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